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Abstract. We apply the hybrid Monte Carlo (HMC) algorithm to the 
financial time sires analysis of the stochastic volatility (SV) model for the 
first time. The HMC algorithm is used for the Markov chain Monte Carlo 
(MCMC) update of volatility variables of the SV model in the Bayesian 
inference. We compute parameters of the SV model from the artificial 
Lj ' financial data and compare the results from the HMC algorithm with 



those from the Metropolis algorithm. We find that the HMC decorrelates 
the volatility variables faster than the Metropolis algorithm. We also 
make an empirical analysis based on the Yen/Dollar exchange rates. 
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>> ' 1 Introduction 



It is well known that financial time series of asset returns shows various inter- 



_ 4- . esting properties which can not be explained from the assumption of that the 

time series obeys the Brownian motion. Those properties are classified as styl- 
ized facts.!.. Some examples of the stylized facts are (i) fat-tailed distribution 

qq ' of return (ii) volatility clustering (iii) slow decay of the autocorrelation time of 

f^) . the absolute returns. The true dynamics behind the stylized facts is not fully 

understood. There are some attempts to construct physical models based on spin 
dynamics [2] -[5] and they are able to capture some of the stylized facts. 

/\i ' The volatility is an important value to measure the risk in finance. The 

volatilities of asset returns change in time and shows the volatility clustering. In 
order to mimic these empirical properties of the volatility Engle advocated the 
autoregressive conditional hetroskedasticity (ARCH) model[5| where the volatil- 
ity variable changes deterministically depending on the past squared value of the 
return. Later the ARCH model is generalized by adding also the past volatility 
dependence to the volatility change. This model is known as the generalized 
ARCH (GARCH) model [TJ. The parameters of the GARCH model applied to 
financial time series are easily determined by the maximum likelihood method. 
The stochastic volatility (SV) model[8 9J is another model which captures the 
properties of the volatility. Contrast to the GARCH model of which the volatility 
change is deterministic, the volatility of the SV model changes stochastically in 
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time. As a result the likelihood function of the SV model is given as a multiple 
integral of the volatility variables. Such an integral in general is not analytically 
calculable and thus the determination of the parameters in the SV model by the 
maximum likelihood method becomes ineffective. 

For the SV model the MCMC method based on the Bayesian approach is 
developed. In the MCMC of the SV model one has to update the parameter 
variables and also the volatility ones. The number of the volatility variables to 
be updated increases with the data size of the time series. Usually the update 
scheme of the volatility variables is based on the local one such as the Metropolis- 
type algorithm [8] . It is however known that when the local update scheme is done 
for the volatility variables which have interactions to their neighbor variables in 
time, the autocorrelation time of sampled volatility variables becomes high and 
thus the local update scheme is not effective. In order to improve the efficiency 
of the local update method the blocked scheme which updates several variables 
at once is also proposed [TO]. 

In this paper we use the HMC algorithmpT] to update the volatility variables. 
There exists an application of the HMC algorithm to the GARCH mo del [12] 
where three GARCH parameters are updated by the HMC scheme. It is more 
interesting to apply the HMC for update of the volatility variables because the 
HMC algorithm is a global update scheme which can update all variables at once. 
To examine the HMC we calculate the autocorrelation function of the volatility 
variables and compare the result with that of the Metropolis algorithm. 

2 Stochastic Volatility Model and its Bayesian inference 

2.1 Stochastic Volatility Model 

The standard version of the SV model [51^j is given by 

Vt = cr t e t = exp(/i t /2)e t , (1) 

h t = n + (j>(h t -i - fi) +T) t , (2) 

where yt — (yi,y2, ■■■■> Vn) represents the time series data, h t is defined by h t — 
In of and 04 is called volatility. The error terms €t and r)t are independent normal 
distributions N(0, 1) and N(0, of) respectively. 

For this model the parameters to be determined are /x, <f> and az. Let us use 
9 as 9 = (/!, <f>, ai). The likelihood function L(9) for the SV model is written as 

/n 
Y[f(e t \ ( x 2 t )f(ht\9)dh 1 dh2...dh n , (3) 



t=i 



where 



yf 



f(e t \at) = (2natr?eM-7^), (4) 

Z(T t 
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f(h t \9) = (2Tra v ) ■ exp( — ). (6) 

2.2 Bayesian inference for the SV model 

In the Bayesian theorem, the probability distributions of the parameters to be 
estimated are given by 

/(%) = ^L(8M9), (7) 

where Z is the normalization constant Z = J L(9)ir(9)d9 and n(9) is a prior disti- 
bution of 9 for which we make a certian assumption. The values of the parameters 
are inferred as the expectation values of 9 given by (9) = J 9f(9\y)d9. In general 
this integral can not be performed analytically. For that case, one can use the 
MCMC method to estimate the expectation values numerically. In the MCMC 
method, we first generate a series of 9 with a probability P(9) — f(9\y). Let 
0W = (0W,6> (2) , ...,6> (fc) ) be values of 9 generated by a MCMC sampling. Then 
using these values the expectation value of 9 is estimated by (9) = \ Yli=i @ ■ 
For the SV model, in addition to 9, volatility variables ht also have to be 
updated since they are integrated out as in eq.©. Let P(9,ht) be the joint 
probability distribution of 8 and ht- Then P(9, h t ) is given by 

P(e,ht)~L(0,h t )n(6), (8) 

where L(9,h t ) = itf =1 f(e t \h t )f(h t \8). 

For the prior n(9) we assume that ir(cr%) ~ (c 2 ) -1 and for others Tr(fi) = 
7r((/)) = constant. The probability distributions for the parameters and the 
volatility variables are derived from eg .((5)) [515] . The probability distributions 
and their update schemes are given in the followings. 

• <7 2 update scheme. 

The probability distribution of er 2 is given by 

PKW^-l^expf-^V (9) 

where A = |{(1 - 2 )(/n - /i) 2 + £"= 2 [^ - M - <f>(ht-i ~ m)] 2 }- 

Since eq.© is an inverse gamma distribution we can easily draw a value 

depending on eq.© by using an appropriate statistical library. 

• [i update scheme. 

The probability distribution of /i is given by 

^M-expj-^Oz-^) 2 }, (10) 

where B = (1 - <j) 2 ) + (n - 1)(1 - <f>) 2 , 

and C = : (1 - 2 )^r + (1 - 0) E? =2 (/*i - 0ftt_i). 

Eq. (fT0| is a Gaussian distribution. Again we can easily update \x. 
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• <j) update scheme. 

The probability distribution of <f> is given by 

P(<f>) ~ (1 - <f> 2 ) 1/2 cxp{-^(0 - |) 2 }, (11) 

where D = -(/h - ^) 2 + I]" =2 (/it-i - A*) 2 and E = Y%=i(ht - ti)(h-i - m )- 

In order to update with eq. (fTTj) . we use the Metropolis-Hastings algorithm |13|14j . 

Let us write eq. (fTTj) as P{4>) ~ Pi(4>)P2{<j)) where 

P 1 (0) = (1-^ 2 ) 1 / 2 , (12) 

P 2 (0)~exp{^(0-|l) 2 }. (13) 

Since P2(</>) is a Gaussian distribution we can easily draw <j> from eq. (|13p . 
Let 4>new be a candidate given from eq. (|13|) . Then in order to obtain the 
correct distribution, <p new is accepted with the following probability Pmh- 



• 



P 



MH 



. \ p{^ new )p 2 {<j)) \ . J rn-^j 

un < — ^—7 ^ ^-f , 1 > = mm < \ / ^- ^- , 1 > . (14 

\P(£)P 2 (^ ero )' / IV (1-0 2 ) ' f l ' 



In addition to the above step we restrict <f> within [— 1, 1] to avoid a negative 

value in the calculation of square root. 

Probability distribution for P(h t ). 

The probability distribution of the volatility variables h t is given by 

P(h!,h 2 ,...,h n )~ (15) 



^P ^ Z^j=lt 2 + 2 e J2cr2/(l-02) Z^i= 



ht- fj.—^ht-i — m)] 2 
2 2ct2 



This probability distribution is not a simple function for drawing values of the 
volatility variables h t . A conventional method is the Metropolis method[13j 
which updates the variables locally. Here we use the HMC algorithm which 
updates the volatility variables globally. 

3 Hybrid Monte Carlo Algorithm 

The HMC algorithm is originally developed for the MCMC of the lattice Quan- 
tum Chromo Dynamics (QCD) calculations pT] where local type update algo- 
rithms are not effective. The notable feature of the HMC algorithm is that it 
updates a number of variables simultaneously. 

Here we briefly describe the HMC algorithm. The HMC algorithm combines 
molecular dynamics (MD) simulations and the Metropolis test. Let f(x) be a 
probability density and 0(x) a function of x = {x%, X2, — ,X n ). We determine the 
expectation value of 0(x) with the probability density f(x) which is given by 

(0{x)}= [o(x)f(x)dx= I ' 0(x)e lnf{x) dx. (16) 
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Now let us introduce momentum variables p = (pi,p2, ■■■,p n ) conjugate to 
the variables x and rewrite eq. (|16[) as 

(0(x)) = i / 0(x)e-? p2+lnfix) dxdp = 4 / 0(x)e- H{p ' x) dxdp. (17) 

where Z — f e~? p dp. H(p,x) is the Hamiltonian defined by H(p,x) = \p 2 — 
lnf{x) where p 2 stands for X)2=iPi- The introduction of p does not change the 
value of (O(x)). 

In the HMC algorithm, new candidates of the variables are drawn by inte- 
grating the Hamilton's equations of motions. The Hamilton's equations of mo- 
tions are solved numerically by doing the MD simulation with a fictitious time. 
To solve the equations we use the standard 2nd order leapfrog integrator. One 
could use improved integrators [IS] or higher order integrators[16 17: if necessary. 

Let (p' , x 1 ) be the new candidates given by the MD simulation. The new 
candidates are accepted with a probability min{l,exp(— AH)} where AH = 
H(p',x') — H(p,x). Since the Hamilton's equation of motions are not solved 
exactly AH deviates from zero. The magnitude of the deviation is tuned by the 
discrete time step size in the MD simulation such that the acceptance of the new 
candidates becomes high. 

For the volatility variables of the SV model, from eq. (fT"5|) the Hamiltonian 
can be defined by 

M* n_f 1 J l f f ''. l e L-t,) [hi-lA 2 y^ [hi - M ~ 4>{hi-i - /i)] s 
H(Pu ht) - ^ - Pi +}^{ Y + ^ }+ 2CT 2 /(1 _^2 ) +2. ^2 

i=l i=l T y Y ' i=2 ^ 

(18) 
where pi is defined as a conjugate momentum to hi. 



4 Numerical Test 

In this section we investigate the HMC algorithm for the SV model with artificial 
financial data. The artificial data is generated with a set of known parameters. 
We try to infer the values of those parameters by the HMC and Metropolis 
algorithms and compares the results. 

Using eq.|T]) with = 0.97,cr^ = 0.05 and fi = — 1 we have generated 2000 
data. To this data we made the Bayesian inference with the HMC and Metropo- 
lis algorithms. The initial parameters are set to </> = 0.5,<7^ = 1.0 and /i = 0. 
The first 10000 samples are discarded as thermalization or burn-in process. Then 
200000 samples are recorded for analysis. The acceptance of the volatility vari- 
ables is tuned to be about 50%. 

FiglD shows the history of the volatility variable /iioo- We use ft-ioo as the 
representative one of the volatility variable. We observe the similar behavior for 
other volatility variables. As seen in FigQ]the correlation of the volatility variable 
from the HMC algorithm is smaller than that from the Metropolis algorithm. To 



Tetsuya Takaishi 




55000 

Monte Carlo history 




55000 
Monte Carlo history 



Fig. 1. Monte Carlo history of HMC (left) and Metropolis (right) in the window 
from 50000 to 60000 Monte Carlo history. 



quantify this we calculated the autocorrelation function (ACF) of the volatility 
variable shown in FigO The ACF is defined as 



ACF{t) 



N 



ELiK?')-<*>)(z(i + *)-<*>) 



(19) 



where < x > and a\ are the average value and the variance of x respectively. 

The autocorrelation time Ti n t of the volatility variables is given in Table 1. 
The values in the parentheses represent the errors estimated by the jackknife 
method. The autocorrelation time is defined by n n t — \ + J2tLi ACF(t). 

The HMC algorithm gives a smaller autocorrelation time than the Metropolis 
algorithm, which means that the HMC algorithm samples the volatility variables 
more effectively than the Metropolis algorithm. 




Fig. 2. Autocorrelation function of the volatility variable ft-ioo for the HMC and 
Metropolis algorithms. We see that the ACF from HMC decreases quickly as the 
Monte Carlo time increases. 
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<i> 


M 


_2 


hum 


true 


0.97 


-1 


0.05 




HMC 

Tint 


0.978(7) 
540(60) 


-0.92(26) 
3(1) 


0.053(12) 
1200(150) 


18(1) 


Metropolis 

'Tint 


0.978(7) 
400(100) 


-0.92(26) 
13(2) 


0.052(11) 
1000(270) 


210(50) 



Table 1. Results estimated by the HMC and Metropolis algorithms. 



HMC 

Tint 



0.960(12) -1.13(8) 0.014(4) 

610(300) 14(6) 1400(800) 55(11) 



Table 2. Results estimated by the HMC to the Yen/Dollar exchange rates. 



The autocorrelation times for the parameters of the SV model are summa- 
rized in Table 1 . The autocorrelation times from the HMC algorithm are similar 
to those of the Metropolis algorithm except for n n t of /i. 

The values of the SV parameters estimated by the HMC and the Metropolis 
algorithms are given in Table 1. The values in the parentheses represent the 
standard deviations of the sampled data. The results from the both algorithms 
well reproduce the true values used for the generation of the artificial financial 
data. Furthermore for each parameter two values obtained by the HMC and the 
Metropolis algorithms agree well. This is not surprising because the same data 
is used for the both calculations by the HMC and Metropolis algorithms. 

5 Empirical Study 

We have also made an empirical study of the SV model by the HMC. The empir- 
ical study is based on daily data of the exchange rates for Japanese Yen and US 
dollar. The sampling period is 1 March 2000 to 29 February 2008, which has 2007 
observations. The exchange rates pi are transformed to r^ = 1001n(pi/pi_i — s) 
where s is the average value of ln(pi/pi_i). The MCMC sampling is performed 
as in the previous section. The first 10000 MC samples are discarded and then 
20000 samples are recoded for the analysis. The estimated values of the param- 
eters are summarized in Table 2. The estimated value of <j) is close to one, which 
means the persistency of the volatility shock. The similar values are obtained in 
the previous studies [ 



6 Summary 

The HMC algorithm is applied for the Bayesian inference of the SV model. It 
is found that the correlations of the volatility variables sampled by the HMC 
algorithm are much reduced. On the other hand we observe no significant im- 
provement on the correlations of the sampled parameters of the SV model. Thus 
it is concluded that the HMC algorithm has a similar efficiency to the Metropolis 
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algorithm and it is an alternative algorithm for the Bayesian inference of the SV 
model. 

If one needs to calculate a certain quantity depending on the volatility vari- 
ables, then the HMC algorithm may serve as a good algorithm which samples 
the volatility variables effectively because the HMC algorithm decorrelates the 
sampled volatility variables faster than the Metropolis algorithm. 
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